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Abstract 

We calculate the annihilation decay widths of spin-singlet heavy quarkonia h c ,hb and rjb into 
light hadrons with both QCD and relativistic corrections at order 0{a s v 2 ) in nonrelativistic QCD. 



With appropriate estimates for the long-distance matrix elements by using the potential model 



and operator evolution method, we find that our predictions of these decay widths are consistent 
with recent experimental measurements. We also find that the 0(a s v 2 ) corrections are small for 66 



states but substantial for cc states. In particular, the negative contribution of 0{a s v 2 ) correction 
to the h c decay can lower the decay width, as compared with previous predictions without the 
0(a s v 2 ) correction, and thus result in better agreement with the measurement. 
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I. INTRODUCTION 



The inclusive annihilation decay of heavy quarkonium is one of the important issues in 
heavy quarkonium physics. Theoretically, it is widely accepted that the heavy quarkonium 
inclusive annihilation decay can be described by nonrelativistic QCD (NRQCD) factoriza- 
tion ffl. In NRQCD, the long-distance effects that cannot be calculated perturbatively are 
described by the long-distance matrix elements (LDMEs), which are classified in the order 
of v, the relative velocity of heavy quarks in quarkonium. Experimentally, a large amount 
of heavy quarkonia data have given more precise measurements for their decay widths and 
branching ratios (for a review and recent results see, e.g.|^]]). Therefore, it is necessary to 
provide more precise theoretical predictions to compare with experimental measurements. 

For charmonium, the inclusive annihilation hadronic decay (into gluons and light quark 
pairs) widths for S, P, .D-wave cc states are all calculated up to 0(a s ) in NRQCD |5|-|TT|. Par- 



ticularly, for the S- wave state r) c , the 0(a s v ) corrections have recently been carried out [12 



where the short-distance coefficients of 0(v 2 ) are calculated perturbatively to next-to-leading 
order (NLO) in a s . After taking the 0(a s v 2 ) corrections into account, the measurements of 
Tj c decay can be described much better in NRQCD [0. For the P-wave state h c , the earlier 
theoretical result at 0(a s ) predicts the hadronic decay width of h c to be about 0.72 MeV0, 
which is a factor of 2 larger than the latest measurements by BESIII, where the central 
value of the total width is about 0.73 MeV and the hadronic decay branching ratio is about 
50% [Q. Thus it is important to do higher order corrections in v to examine whether the gap 
between theoretical predictions and experimental measurements can be reduced, and this is 
also an interesting test of the validity of NRQCD factorization. 

For bottomonium, the bb system, the value of v 2 is about 0.1, which is much smaller than 
v 2 0.3 for charmonium. It is then expected that the v 2 expansion should be better for 
bottomonium, thus the v 2 correction for bottomonium is more solid to check NRQCD fac- 
torization. Recently, the process hb(lP) —> rjb(lS)j is measured by the Belle Collaboration. 
They find the rjt, decay width to be about 12.4 MeV and the decay branching fraction of 
B[h b (lP) r}b(lS)i] = 49.2 ± 5.7+|3%[|]. It is tempting to try to explain these data in 
NRQCD. 

In this paper, we will perform the 0(a s v 2 ) calculations for the spin-singlet P-wave char- 
monium h c and bottomonium hb, and also for the spin-singlet S'-wave bottomonium rjb- We 
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find these corrections are important to understand the measured data. The rest of this 
paper is organized as follows. In Sec. [TT] we briefly introduce the NRQCD factorization for- 
mulism in heavy quarkonium annihilation decays. Then we describe some technical method 
in calculating 0(a s v 2 ) short-distance coefficients in Sec[m]. The results for S-wave and P 
wave-states including real and virtual contributions are presented in Sec. |1V With these 
results and appropriate estimates of the LDMEs, we discuss the related phenomenology in 
Sec. M. Finally, we give a brief summary in Sec. N% 



II. NRQCD FACTORIZATION FOR QUARKONIUM DECAY 

In this section, we introduce the NRQCD factorization formula for the rates of spin- 
singlet heavy quarkonium {j] c ^ and h Cy b) decays to light hadrons. The inclusive annihilation 
decay width of heavy quarkonium can be factorized by the following formula |I| 

nH) = J2 2Imf ^ A) (H\OM\H), (1) 
rrin 

n V 

where Im/„(^) is the short- distance (SD) coefficient which can be perturbatively calculated 
with full QCD Lagrangian. The long-distance matrix elements(LDMEs) (H\O n (^\)\H) 
involve non-perturbative effects and are classified by the relative velocity v of Q and Q. 

The NRQCD Lagrangian can be derived by integrating out the degrees of freedom of 
order tuq, the mass of the heavy quark 



-NRQCD 



Aight + -Aeavy + SC. (2) 



The heavy part of the Lagrangian describes the motions of (anti-)heavy quark in spacetime 
and is given by 

D 2 D 2 

Aieavy = ^( iD t + 7, )V> + X^ (l'A + „ )X (3) 

zrriQ iniq 

where ip(x) denotes the Pauli spinor field that annihilates (creates) a heavy (anti-) quark, 
and D t (D) is the time(space) component of the gauge- covariant derivative D^. The light 
piece of the Lagrangian reads 

Aight = -~Tr G^G, V + ( 4 ) 
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where G^ u is the gluon field strength tensor, q is the Dirac spinor field of light quarks and 
71/ is the number of light flavors. The bilinear Lagrangian term which contains the order v 2 
corrections can also be reproduced 

5L bilincar = _^t( D 2)2^ + ^^t( D .0 E _0 E . D ty, 

8m Q 8m Q 

+ -%-ibUiT)xg~E- oExzD)-cx^ + -^-^(gB-a)ib 

8m 2 Q • 2m Q 

+ charge conjugate terms, (5) 

where E i = G 0i and B i = \ e ijk & k are the electric and magnetic components of the gluon 
field strength tensor G^ u , and q = 1 + 0(a s ), i = 1, 2, 3, 4 are the dimensionless coefficients 
corresponding to each operator. 

In order to describe the annihilation decay of quarkonium, a set of local 4-fermion op- 
erators Oi which appear in Eq. ([[]) are generated. For example, the operator ^XXT can 
annihilate a QQ pair in the configuration. In our case, for the 0(a s v 2 ) calculation 



of spin-singlet quarkonium decay, the power counting rules [p~3|Hl5|| give the following six 
operators and LDMEs in Eq. (H): 



for S wave, 



for P wave 



and 



OfS®) = ^XXH, (6a) 
VCS l o ] ) = ^ t xx t (-^)¥ + h.c, (6b) 



(7a) 
(7b) 

(7c) 

(7d) 



(0{ 2S+l L [ j' 8] )) H = (H\0( 2S+1 L [ 1' 8] )\H), (8a) 
{v( 2S+i l M ))h = {H \V( 2S+1 L [ j' s] )\H). (8b) 













^ + h.c, 


oCp1 1] ) 


g 

= ^ f ( — ^)x-x\- 


-% 




•o- 

= ^ ( — ^)x-r 


(- l —) 3 ?p + h.c. 
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By counting the dimension of the above operators JT|, we can explicitly rewrite factorization 
Eq. ([I]) for and x Pi states, 

F(H(%) -+ LH) = ^l(OCsS ] )hs + ^^(VCS^hso, (9a) 

m Q m Q 
m Q m Q 

+ ^^-(OCP?))^ + ^p-p^))*. Ob) 



Through the above factorization formula, the work is simplified by matching full QCD 
with NRQCD to get the imaginary part of the short- distance (SD) coefficients F and G 
perturbatively. The skeleton of the matching procedure is given by 



lmM(QQ^QQ) = T 2 lm { n{ Jl k \ QQ\On^A)\QQ) , (10) 

pert QCD ^ m7? NRQCD 
n V 

and the procedure will be discussed in detail in next section. 



III. DETAILS IN FULL QCD CALCULATION 
A. Kinematics 

We investigate the 4-momentum of heavy quarkonium in its rest frame. It is customary 
to decompose the momenta of Q and Q in the following form, 

PQ = \ P + ^ ( lla ) 

pq = \p-q, (lib) 

where 



P = (2£ q ,0), (12a) 
q = (0,q), (12b) 



where _E q = \ uOq + q 2 and satisfies the relation P-q = 0. 
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The treatment of final state phase space of 0(a s v 2 ) corrections is slightly different from 
ordinary calculations (i.e. leading order of v calculation). We make the the following rescal- 



ing transformation for all external momenta in amplitudes |T3, [L6 



P ->• P'^-, (13a) 
m Q 

k f ->• k'^-, (13b) 

but maintain the relative momentum q and loop integral momentum I unchanged. Once we 
take such a trick, the q 2 dependence both in phase space and current factor (i.e. 1/(2M) 
where M is the quarkonium mass) can be absorbed into amplitudes, then we can safely take 
q — > in these terms and only expand q in amplitudes. Here we should note that this trick 
can only work in the case that all the final state partons are massless (i.e. gluons and light 
quarks), because for the massive partons the on-shell relation under rescaling do not hold. 



B. Covariant Projection Method in D-Dimension 

We use an equivalent but more efficient method, i.e. the covariant projection method, 
to calculate the imaginary part of SD coefficients in Eq. ( pa]) and fl9bf) . In order to get 
spin-singlet color-singlet and color-octet QQ decay amplitudes, we take the following spin 
and color projectors onto QQ quark lines |fl7f1 : 



n n 



2V2( Eci + m Q ) { 2 + * + mQ) 8^ ( 2 



m Q ), (14) 



and 



Ci 



V2T C . 



(15a) 
(15b) 



Since the projected amplitudes also depend on q, we should make Taylor expansion in powers 
of q to the required order, 

dM{q) 



M(q) = M(0) + 
+ 



dq' 

1 d 3 M(q) 



3! dq a dq?dq^ 



9 =o 



1 d 2 M(q) 
q =o q 2! dq a dqP 

q a q^q" + • • • , 



g=0 



(16) 
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and then make the replacement: 



q a qp -> 



D-l 



-n 



a/5, 



1 I 

q q- q 

QaQ^Qx -> ^ - 1 (n a/3 n 7A + n a7 n^ A + n aA n 7/3 ), 



where 



n 



a/3 



-9aB + 



pi pi 

4m 2 Q ' 



(17a) 
(17b) 

(18) 



and q' is the relative momentum corresponding to the complex conjugate of the amplitudes 
Ai'(q') and P' is the rescaling heavy quarkonium momentum. For example, the third 
derivative term of M. times the first derivative term of Jvft gives the squared amplitudes 
term 



1 d 3 M(q) dM\q' 



=o dq 



i\ 



q'=0 



q a q fi q 1 4 x 



3! dq a dq?dq^ 

1 qV^ m n + n n + n n ^ d " M ^ 



dM^{q') 



=0 dq 



IX 



q'=0 



(19) 



and this term contributes to the SD coefficient of Gi^Pl ) in Eq. (BE). 



IV. PERTURB ATIVE QCD RESULTS OF SHORT-DISTANCE COEFFICIENTS 

Our techniques of calculating the SD coefficients at 0(a s v 2 ) are as follows. First we 
generate Feynman diagrams and amplitudes by FenyArts [ffSl [H5J, and then calculate the 
squared amplitudes by self-written Mathematica codes, and the phase space integrals are 
calculated analytically by the method presented in Ref. ||. Ultra-violet (UV) and infra- 
red(IR) divergences are regularized by dimensional regularization. The renormalization of 
UV divergences is taken by re-defining the renormalized heavy quark mass itiq, heavy quark 
field ipQ, light quark field ip q and gluon field in the on-mass-shell scheme(OS), and the 
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FIG. 1: Born level Feynman diagrams for 
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QCD coupling constant g s in the MS scheme, 
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47T 
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5Z 



MS 



' 2 4tt e 



-2CU) 
1 



(20a) 
(20b) 
(20c) 
(20d) 
(20e) 



where N e (rriQ) = (^^) e r(l + e) is an overall factor, ft = yC^ — ^T F rif is the one-loop coef- 
ficient in the ft function of QCD, n/ is the active quark flavors, and // r is the renormalization 
scale. 



A. Short-Distance Coefficients of S-Wave Quarkonium in Hadronic Decay 

Our leading order calculations in a s give the Born level decay width and its relativistic 
correction respectively as 

rBomC^ ^) = 3(47T« S ) 2 ^$( 2 )(1 - e)(l - 2e) V ° , (21a) 

r£ n (^ ] ^) = "^^^W^ 99), (2ib) 

where $( 2 ) = ^(jjsY r(2-2% * s ^ e total two-body phase space in D dimension and M = 
2m,Q y/l + -^r is the quarkonium mass including the relativistic correction. The two Born 
diagrams are illustrated in Fig. [I]. 

The next-to-leading order calculations include real and virtual corrections. For S-wave 
Fock states (i.e. x Sq and 1 S^), UV divergences will be canceled by counter-term diagrams, 
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FIG. 2: Virtual correction Feynman diagrams for iff? 1 
state. 



1 f8l 

55. Only Sq contributes for P-wave 



and IR divergences will be canceled each other by the two corrections in the leading term in 
v 2 , but some soft divergences remain in the relativistic corrections. The cancelation of such 
additional divergences will be presented in next section by calculating NRQCD LDMEs at 
1-loop level. The contribution of virtual plus counter-term corrections is 

rVirtual^i 11 99) =— r B om( 1 ^ 1 99)fe(m Q ){[-^ 

19vr 2 - 44)] 



e 



1 4m 2 n 

+ 36 ( - 6/3 ° ln( 7f 
^ 2 r 4 1 



_<r_ 

m 2 [ 3e 2 



4n 



27 



971 

e 



(22) 



1 



324 



-72/3 ln(— ^)+8n/ 



267tt- 



- 280)]}, 



2 

where / £ (mq) = (^r)T(l + e). Some selected Feynman diagrams are shown in Fig. 0. 

Q 

The real correction contains two sets, where the final states of one set is three gluons and 
the other one is qqg. The diagrams are shown in Fig. ^| and Fig. ^ and the contributions to 
decay width are shown respectively as 

TCS [ o ] 999) =— rBornC 1 ^ 1 ^)/e(mo){[i + ^- + ^(724 - 69vr 2 )] 
7r l e z 6 e 72 

4 1 3 437-42tt 2 
3?~7 



+ 



m. 



r^ 1 -+ gqg) =^-r Bom ( 1 4 1] -+ 99) 



27 



(23a) 



fe{m Q ) 



2 7T 



r i + e r 1 



2 1 16 q 2 8 1 86 



(23b) 
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FIG. 3: Real correction Feynman diagrams for ^ ggg. The crossed diagrams have been 
suppressed. It is the same for both l pj^ and Fock states. 







— 1575757515^ 
v 

FIG. 4: Real correction Feynman diagrams for ^j, 1 ' — > qqg. Only ^q 8 ' contributes for the P-wave 
state. 

Combining Eqs. (|21"D, ( f2"2"|) and ([23|), we obtain the hadronic decay width with both QCD 
radiative and relativistic corrections at NLO of ^So heavy quarkonium, 

Tqcd^o -»- LH) =r Bom ( 1 4 11 -»> 99){ [I + ^/e(m o )l(-36/3 ln( 4m "' 



2 



7T " ' ^72' V fJh 

64n f - 937T 2 + 1908)] - \\ [l + —f e {m Q ){-\- (24) 



1 477Z 2 

+ ^(-72^o H—r) ~ 164 ™/ - 237 ^ + 4964))] }. 



1 , , , 4m< 

We note that our results agree with the previous work for 0(a s v 2 ) correction [[T^j and 0(a s ) 
correction || |TT|. By correctly comparing our results with Ref. |12|| , a slight difference of 
two body phase space $2 between them can be found. In Ref. |T2| $ 2 is defined by removing 
the q 2 dependence into the coefficients, so our individual virtual and real parts Eq. (^) and 
Eq. ( p3"D look different from the results in Ref. JL2] but essentially they are equivalent. The 
total NLO result Eq. (|24[) is explicitly the same, independent of the definition of $2- The 
correct repetition of the hadronic decay SD coefficients of ^So heavy quarkonium enables our 
extended discussion from charm quark system to bottom quark system (i.e. ■%) and also 
partly checks our codes when dealing with P-wave heavy quarkonium. 
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B. Short-Distance Coefficients of P-Wave Quarkonium In Hadronic Decay 



The procedure in calculating the hadronic decay width SD coefficients of x Pi heavy quarko- 
nium is similar to but more complicated than the ^So case. Additional simplification can be 
taken by imposing C (charge) parity conservation of QCD in selecting Feynman diagrams. 
A straightforward result is that C parity conservation prohibits Fock state, which has 
C = —1, to decay to two gluons, whose C = +1, no matter they are real or virtual. The 
results are shown as follows. 
At the Born level, 

r B orn( 1 4 81 99) = ^(47ra s ) 2 ^|$ (2) (l - e)(l - 2e)(OCS [ * ] ))?°™ (25a) 

12 m?i 



AcP] , „„\ _ 2 ( 2 - e ) V 



2 



r BomC S o l -+ 99) = -^r^BomCS 1 ^ -X 99), (25b) 
" Ulli 6 — 2e wig 

For NLO corrections, 



3a* „ ,i „rsi , „ , s f\ 1 n/ - 21 1 

9 e 



rVirtuaK^O 81 -> =^ £r Born( 1 4 81 ^ ^)/e(m Q )|[-- 

1 4m 2 
+ -(-12/3 ln(^) + 297r 2 -16)] 
72 ^ 

q 2 4 1 An f - 115 1 

+ 27 7 

1 4t?2 2 "I 

144/3 la( — + 16n/ + 345tt 2 - 992)] }, 



(26) 



628 v //; 

r( 1 4 81 ^) =^r Born ( 1 4 81 -+ ^)/ e (mg){[i + 7 -- - vr 2 + ^] 

2 4 1 4 554 -45tt 2 ) 
+ ^~3? "7 27 J /' 



l(^o 2 /Borneo ^^r(l + e)r(l-e) [ 36 9 + m| ( 9 6 + 27 )J ' 

(28) 

40a? . , w _ , f r 1 7tt 2 5, 



T( 1 P 1 [1] ^) =^/ e (m )(87r$ 2 ){[-i + 



24 3 

q* 291 4216 - 555tt 2 -| (C^if 1 ))]*"" 1 (29) 
+ m^l5e + 900 J J 2A" c m^ ' 
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Summing over the above results we get the total hadronic decay width, 

r QC D(^i LH) =r Born ( 1 ^ 1 gg)[[l + ^ /e ( mQ )(-I/3 Q ln^ 



t4 



9 ' 24 J 3m 2 Q l 7T JeV gn 12 e 
+ _(-l44A,ln(— r) - 328n / - 7357r2 + 12304 ))] 1 ( 3 °) 

ZOO [I* J 

40a? . , w . T j r 1 7tt 2 5, 



+ ^/e(m Q )(8vr$ 2 ){[ 

q* 291 4216 - 5557T 2 -j (ggfj))^ 



m^ L 15e 900 J J 2iy c m^ 

C. Evaluating NRQCD LDMEs And Matching Full QCD Results 

In Eqs. (p4|) and (|30"D there exists explicit IR divergences. This is because, in principle we 
should evaluate both SD coefficients and LDMEs at the loop level and then take the matching 



procedure (10). So all the Born LDMEs appearing in Eq. fl24|) and ( p0| ) should be replaced 
by one-loop LDMEs by calculating the loop-corrections for them. Once we make such 
a replacement, all IR divergences should be canceled and the results become infra-red safe 
quantities. The self-energy contributions which connect Born LDMEs to their corresponding 
relativistic ones are first calculated in the Appendix of Ref . |IJ . The intersecting diagrams 
which describe the El transition between and states at 0(a s v 2 ) in this work is 
new. The detailed calculation is presented in Appendix Here we give the relevant results 



in dimensional regularization and renormalize them in MS scheme 

,2 1 .,2 



(OesV))?°'° -><OCS™))i£>{l - \^^Um Q )[\ - K£|)]}, (31a) 



16a s „ . , r r l , . iA 



5m 2 [ e l 4m2 jJ J 2N,m 2 : 



5m 2 , e 4m 2 , ) 2N c m 2 Q 
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where /ia is the factorization scale. Substituting them into Eq. (p4|) and (p0[), and consid- 
ering the relation 

(res^r = ow^S" 1 , ( 32a ) 

(P(^ ] ))P P r = q 2 (0(^ ] ))g? m , (32b) 

^ipJl^Born = ^^(XpW^Born (32c) 

we get the SD coefficients for heavy quarkonium hadronic decay of S-wave and P-wave states 
by matching full QCD and NRQCD, 

Anal 



ictib 



FfS [ * ] ) 
G( 1 4 81 ) 

G( 1 P 1 [1] ) 
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a a 1 4m?> 
1 - — =r 36/3 ln — ^ + 64 % + 93tt 2 

7T 72 /i* 



1908) 



AAna 2 s 
3 9 



1 



7T 144 



192 ln( 



/4 

4m Q 



+ 72/3 ln( 



4m|, 



+ 164n/ + 2377r 2 -4964]} 5 



57ra 2 

' — 6~ 



1 



1 

— — (36/3 ln(—^) + 64n f + 129vr 2 - 244£ 

7T 72 



A5na 2 s 
3 6 



1 - 



7T 288 

+ 32871/ + 735tt 2 - 12304] 
5a? 



' [168 ln(-^ 



"4m^ 



+ 144/3 ln(- 



4m 2 



/'r 



486 

„3 



7(7r 2 -16)-241n( 



4m 2 / 



3645 



1740 ln( 



j4_ 

4m 2 - 



- 555tt 2 + 9236 



(33a) 

(33b) 

(33c) 

(33d) 

(33e) 
(33f) 



where F's and G's are defined in Eq. fl9a|) and (9b"). 



The SD coefficients of agree with those in Refs. |, g |TI|, @ 0, that of ^ 8J and 
at leading order in v 2 are also agree with previous results in Ref. ||. The relativistic 
corrections G^S^) and G( 1 P 1 ) are primarily new results in this work. Based on these 
results, we will further analyze the decay of 1 Sq and l P\ heavy quarkonium into light hadrons. 



V. PHENOMENOLOGICAL DISCUSSIONS 



A. Estimating NRQCD LDMEs 



As shown in Eqs. ( |9aD and (9B) in previous section, the decay widths of heavy quarkonium 
are split into two parts, one is the SD coefficients which can be perturbatively calculated 
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in the framework of NRQCD factorization and the results are listed in Eq. fl33|), the other 
is LDMEs containing non-perturbative information and cannot be calculated precisely at 
present. For ^Sq quarkonium there are two LDMEs and for x Pi are four. In Ref. |T!I the 



LDMEs of T] c charmonium are determined by using the Cornell potential [|2D| and one of the 
experimental results of r LH (?7 c ) and r 77 (r] c )[^T[], and then predict the other. In the present 
work, since there are no enough experimental information to determine all involved LDMEs, 
we will estimate the values of (i-e. rjb) and (i.e. h c and hb) LDMEs separately. 

For rjb quarkonium, the situation is similar to Ref. [12], but lacking the experiment input 
of the decay width to two photons r 77 (^). In this case we may first determine (OQSq )) m 
from potential model. Here we use the Buchmuller-Tye(B-T) potential model and 



Cornell(Corn) potential model [^Oj results as input. They are p2l p3| 



(Offlf = ^\Rf T (0)\ 2 = 3.093 GeV 3 , (34a) 
(OCS^r = (0(*>?S = 3.071*3 GeV 3 . (34b) 

Here we use the heavy quarkonium spin symmetry, which holds in leading order of v, and 
higher order corrections for spin symmetry are expected to be small and less than the 
uncertainty of potential model estimates themselves, to relate rjf, and T(IS'). In order to 
determine (V( l S^)) Vb , we define JT2], [TO 



/ 2\ _ ( V ( S o Thy, /oc\ 

although (v 2 ) Vb can not be understood as the expectation value of v 2 in potential model, it 
can be estimated from the Gremm-Kapustin relation: 

^2^G-K _ m Vb ~ 2 m pole ^gQ-j 

Here we choose m po i e = 4.6 GeV for b quark and m % = 9401 MeV from recent experimental 
measurement. Then we get {v 2 ) Vb = 0.044, which is close to the potential model estimate 
v 2 ~ 0.05 - 0.1. 

For h c quarkonium, the estimate is more complicated since we need to determine four 
LDMEs. (O^Pi )}h c is determined by the B-T potential model and 

(VCPi 1] )) hc = (v 2 ) h m 2 c (OCPi 1] )) hc » (v%ml(OCP?)) hc , (37) 
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'Ic 



where {v 2 ) ric = 0.228 from Ref. [0. Here we have tentatively assumed {v 2 )h c ~ (v 2 ) 
The remaining two color-octet LDMEs are determined by the operator evolution method 
(OEM) J?§ [25J. From Eq. ( |31b|) we get the evolution equations 



772 ~~ ~777 772 h 



dfi\ 9tt m 2 Q 9vr 2N c m 2 Q 15vr 2iV c mJ 

2 d(VeS [ * ] )) 16a s {V?P?)) 



(38) 



d/4 ' 971 2N c m Q 

by integrating them over the range from n\ = m c v ~ 0.8 ± 0.2 GeV to 2m c for the typical 



charmonium NRQCD scales (here we neglect the initial values at m c v [P5|). Here we set m c 
to be its pole mass 1.5 ± 0.1 GeV, and the results are 

C C 

<Wi W )>a. - = 0.8061°;^ MeV, 

{Wst ] )) hc = { { ° 2 ))hc = 2.5Utlf 2 t MeV, 

(p ( i 5 M )>hB ^ gCgglk = 0.479±SiSS MeV. 

The above values meet the NRQCD velocity scaling rules and it is appropriate to determine 
CO matrix elements by OEM. The uncertainties are from variations of m c , (v 2 ), \R'p{0)\h c for 
color-singlet (CS) matrix elements and also from the uncertainty of fi\ for color-Octet (CO) 
matrix elements. The related phenomenology about the uncertainties that might induce the 
variation of the decay width will be discussed in detail below. 

Using the same method we can determine the LDMEs for hj,, they are 

(0( 1 P 1 [1] )), 6 = 755.5iljkeV, (P( 1 P 1 [1] )), 6 = 33.2^?keV, 
____ _ (40) 

(OCS [ * ] )) hb = 441.6^keV, (VCS^)) hb = 18.7±^keV. 
Here we choose = 4.6 ±0.1 GeV, NRQCD scale fi\ = ra^v ~ 1.5 ± 0.2 GeV and set 
{v 2 )h b ~ ( v2 )r]b f° r ^he same assumption as h c . With these results of LDMEs, we can get 
the theoretical decay widths and compare them with the experimental data. 

B. T(r] b -> LH) 

With the LDMEs results given above (here the values (O^S^)}^ and (v 2 )^), we can 
estimate the hadronic decay width of rjt- First we fix both the renormalization scale [i r and 
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factorization scale //a to be 2m& and consider the difference choice of LDMEs. At this scale 
the decay width can be rewritten as 

T( Vb LH) = 441.1j*g x lO- 3 ^ 1 ^ 1 )),, - 664.8j*g x lO^CS? 1 )),, (41) 



where 



I C J|,\ _ (^(^Q 1 ))*ft 

% " 2iV c mg ' 



c b 



(42) 



The LDME (C( 1 ^ % 6 is determined by Eq. ( piaD , (0(% > ))^ L = 24.361}$ MeV and 
(0( l S [ o ] ))%°™ = 24.18^ MeV. (P^U is smaller than (0(^ ] )) m and get 

(WS [ o ] ))n T = l-07±J5g MeV and (n 1 ^))?™ = 1.06j£gg MeV. 

These combinations of LDMEs give the hadronic decay width of T]b, 

T(r] b -»■ LH) B " T = 10.03+°:45 MeV ' ( 43a ) 
r(^^LH) Corn = 9.96±g;|l MeV, (43b) 

and the /i r dependence is plotted in Fig. ||. Here we see that the relativistic corrections 
only contribute about 10% to the total decay width and may not be important. 

In order to compare with the experimental data of rjb total decay width ||, some ap- 
proximation should be mentioned. The difference between the total decay width and the 
hadronic one of rjb in QCD is mainly from the process r]b — > ccX, whose Feynman diagram 
are given by just replacing the final states qq pair in Fig. [| by cc pair. Since it has no tree 
level contribution and is suppressed by the final state phase space, we can safely estimate 
the decay width of this process to be of the order a s (l — ^) 2 ~ 5% less than the decay 
width into light hadrons. This value is within our error band and can omitted. So we have 
T(r]b — > LH) w T tota ^(r]b) ~ 9 — 12.5 MeV. This value is consistent with the experimental 
data T ex P(r] b ) = 10.8±|?±^ MeV §. 



C. T(h c LH) 

The numerical values of SD coefficients for hadronic decay width of h c are 
T(h c LH) =328.7™ x lO^fS®))^ - 39.6l 3 j x 10^ (0( l Pl 1] )) hc 

(44) 

- 446.0™ x 10- 3 (^(^ ] )}, c + 92.4^ x lQ-^P?))^ 
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MGeV) 



FIG. 5: fi r dependence of T(j]b — > LH). LO represents values without QCD and relativistic 
corrections, NLO* includes QCD corrections but only at leading order in v, and NLO takes into 
account all contributions at 0(a s v 2 ) including both QCD and relativistic corrections at next-to- 
leading order. The LDMEs are taken from the B-T potential model and the Gremm-Kapustin 
relation. Here fi\ is set to 2m c . 

where we also set /i r = fi\ = 2m c to eliminate logarithms and define the modified LDMEs 
as in Eq. ([39|). From these values we can see that the color-octet part, which contains larger 
SD coefficients, is the primary contribution to the decay width. We also note here that for 
cc system the relative velocity squared v 2 ~ 0.3 is not small so the v 2 power counting law 
does not work very well to classify the LDMEs. In fact we can see from Eq. (^) that the 
relativistic corrections are substantial as compared with the leading order contribution. The 
decay width including the total and each LDME contributions are listed in Table [p. 



TABLE I: T{h c — > LH) expressed with the contributions of each LDME. 





(oes [ * ] )) hc 


{Oi}P?)) hc 






Total 


T( 2S+1 l)j — >■ LH)(MeV) 


0.83^1 


_n 14+0-03 


n 91 +o.09 

U - Zi -0.14 


07+ - 04 

u - u ' -0.03 


u - oo -0.19 



The results show that the relativistic correction is important for h c hadronic decay, es- 
pecially the color-octet correction. We can also see the contributions from different orders 
of a s and v, which are listed in Table O. The /i r dependence are also shown in Fig. [| 
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TABLE II: T{h c — > LH) expressed with contributions at various orders of a s and v. 





a y 


alv° 


a y 


ay 


Total 


r(/i c ^LH)(MeV) 


o 4o +0 - 18 




-0 10+ 004 
u - iu -0.07 


-0 04+ - 03 


u - oo -0.19 




FIG. 6: n r and z/a dependence of T(h c — > LH). The upper plots are for \i T and lower ones for fj,\. 
From left to right the plots are shown for LO, NLO* and NLO respectively, where NLO* includes 
0{a s ) but excludes 0(a s v 2 ) corrections. 

Here we can see that the contributions of the first two terms to the decay width at 
leading order in v are both positive, whereas the relativistic corrections are negative. These 
relativistic corrections play an important role and pull down the theoretical values of h c decay 
width by partially canceling the contribution at leading order in v. In order to compare with 
the experiment data [|], we should evaluate the El transition decay width T(h c — > rj c + 7) 
also up to v 2 . Ref. estimated the transition decay widths but only at leading order in v 2 
by using the spin-symmetry relation between the spin-singlet and triplet P wave charmonia 
(weighted with a factor of 2 J + 1), 

T(h c ^ 1Vc ) = — — ^(2J+1) ^jy 3 . (45) 

And the obtained El width is 615 ± 29 keV with the spin-triplet input in Eq. ([45]) extracted 



from PDG Data [21]. This result is consistent with the theoretical calculations at leading 
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order in v [[42]. However, if we consider the v 2 corrections, the spin-symmetry relation will 
not hold at this order any more, and instead we should consider the relativistic corrections 



for wave functions. Ref. showed that the width of h c — > ^rj c can be reduced from 650 KeV 
to 385 KeV by relativistic effects. Subsequent studies in various potential models [l|3|-05 



also observed these relativistic effects on this El transition width, ranging over 354-323 KeV. 
In this paper we choose the value T(h c — > 7^ c ) = 385 keV from Ref. |X2 . 



These two decay channels (to LH and to 7^ c ) are the dominant ones of h c . The summation 
of the two contributions gives the theoretical total decay width T th (h c ) = 0.93i£ig MeV and 
the branching ratio B^(h c — > r\ c + 7) = 4llg°%. The total decay width falls in the error 
band of experimental data r exp (/i c ) = 0.73^Q28 MeV, while the branching ratio is a little 
lower than the experiment £> exp (/i c — > rj c + 7) = 54.3 ± 6.7 ± 5.2%. However, if we ignore 
the relativistic corrections to the hadronic decay width, the total width will increase to 1.07 
MeV and the El transition branching ratio will be pulled down to 36%. Therefore, it is 
evident that the relativistic corrections play an important role in the h c decay and can lead 
to a better agreement between theoretical estimate and the experimental data. 

The uncertainties of hadronic decay width are mainly associated with the determination 
of color-octet (CO) matrix elements (see the second column in Table. ||). In OEM they are 
evaluated by integrating out the leading logarithms from /ia to /iA and are then proportional 
to ln(— ft * s known that when the scale is close to Aqcd, the error of a s becomes 
large. Therefore, the CO matrix elements are sensitive to the value of q; s (//a ). In Fig. [7] we 
illustrate the /ia dependence of hadronic decay width in comparison with /i r , /ia dependence 
in Fig. ||. We see that the decay width is indeed sensitive to the value of /ia when [i\ < 1.5 
GeV. So we need to invoke other processes to determine it. If we assume that the initial 
scale fix is universal in all P-wave charmonia decays, we may determine it by e.g. Xcj 
LH(77)(J = 0, 2), whose decay widths have been precisely measured, and this kind of studies 
may be left as our future work. The other underlying uncertainty is the color-singlet matrix 
element ((9( 1 P| 1 '))/ lc ~ 3.5 MeV which is determined by potential model in our calculations. 
We note that in Ref. || it was extracted from the data of Xcj decays and the value is 2.30 
MeV, but this was the leading order result in v, where only two LDMEs are involved. At 
order v 2 many more LDMEs appear and need to be determined, which can not be extracted 
from limited experimental data at present. Nevertheless, if we take this smaller value 2.30 
MeV as our input and keep relations Eqs. fl3"T| ) and ( |3lf ) unchanged, the hadronic decay 
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LO NLO* NLO 




0.6 0.8 1.0 1.2 1.4 0.6 0.8 1.0 1.2 1.4 0.6 0.8 1.0 1.2 1.4 

fA,(OeV) ^ A „(GeV) ^a„(GcV) 



FIG. 7: /xa dependence of T(h c — > LH). From left to right the plots represent LO, NLO* and 
NLO respectively, where NLO* includes 0(a s ) but excludes 0(a s v 2 ) corrections. Here /i r and /ia 
both are set to be 2m c . 

width will be pulled down to 0.36^o^ MeV. Then the total decay width and El transition 
branching ratio will become 0.741q^ MeV and 52^j%, which fit the experiment central 
value very well. Despite of all these uncertainties, we see that the negative contribution of 
relativistic corrections is certainly in favor of explaining the small decay width of h c . 



D. T(h b LH) 

The analysis is by analogy with h c decays except that the relativistic corrections are not 
so large for bb system as discussed in r\ h decay. The SD coefficients and the decay width are 

T(h b -> LH) =(150.2 ± 2.1) x lO- 3 ^ 1 ^ 81 ))^ - (15.3 ± 0.3) x lO- 3 ^^ 11 )}^ 

_ _ ( 46 ) 

- (203.3 ± 2.9) x lO- 3 ^ 1 ^ 81 ))^ + (35.8 ± 0.6) x 10- 3 (P( 1 P 1 [1] )) hi) , 



TABLE III: T(h b — > LH) expressed with contributions of each LDME. 





(oes [ * ] )) hb 


(oepi 1] )) hb 






Total 


Y{ 2S + x L { f -> LH)(keV) 


66.33±l;fg 


-11.581SS 


q 01+O.6I 
" 3 -° 1 -0.72 


i.i9ig:g 


i q+8.62 



TABLE IV: T(h b — > LH) expressed with various orders of a s and v. 





a y 


a y 


a y 


a\v 2 


Total 


T(h b ^LH)(keV) 




1 7 4Q+ 3 - 65 
11 - 4y -3.04 


o i 1+0.34 


U - Oi -0.21 


co 1 q+8.62 
<JZ.±0_ 7 2g 
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FIG. 8: fi r dependence of T(h b —> LH). From left to right the three plots represent LO, NLO* 
and NLO respectively, where NLO* includes 0(a s ) but excludes 0(a s v 2 ) corrections. 



And the fi r dependence is plotted in Fig. [|. From Table IV we see that the v 2 corrections 
are much smaller than the total decay width. The El transition decay width for h b is 



evaluated in the NR |2q| , GI ||4| and Screened-potential models P6[ . The results are listed 
in Table 0. 

TABLE V: T(h b -)• r] b + 7) and B{h b -)■ rj b + 7) in NR, GI and Screened-potential models(SNR 
is calculated using the zeroth-order wave functions while SNRi using the first-order relativistically 
corrected wave functions) 





NR 


GI 


SNR 


SNRi 


T(h b -> rjb + 7) (keV) 


41.8 


37.0 


55.8 


36.3 


rtotal(^) ( keV ) 


93.9 


89.1 


107.9 


88.4 


B{h b r] b + 7) 


44.5% 


41.5% 


51.7% 


41.0% 



Compared with the new experimental data B exp (h b (lP) -> r] b (lS)^) = 49.2=L5.7±|;|% [|, 
the NR and SNRq results fit it well while the GI and SNRi ones are somewhat lower. 



VI. SUMMARY 

We have calculated order a s v 2 corrections for the annihilation hadronic decay widths of 
spin-singlet heavy quarkonia 7] b , h c and h b within the framework of NRQCD. The short- 
distance (SD) coefficients are calculated by the method with matching procedure, and the 
long-distance matrix elements (LDMEs) are estimated by using the potential model and 
operator evolution methods. For the h c decay, we find that 0{v 2 ) and 0(a s v 2 ) corrections 
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contribute large but negative values to the decay width, which substantially reduce the decay 
width calculated in the leading order in v 2 . It shows that relativistic corrections play an 
important role in hadronic decays of cc system, and can improve the theoretical results as 
compared with experimental data. Our calculated total decay width is V (h c ) = 0.931q 'H 
MeV and the branching ratio B^(h c — > r\ c + 7) = 41^g°% with the color-singlet matrix 
elements calculated by potential model, and T til (h c ) = 0.74+QJ5 MeV and B^Qic — > r\ c + 
7) = 52^% with the color-singlet matrix elements extracted from data (but at leading 
order in For rjt and hb decays, we have calculated their hadronic decay widths and 

found that T(rj b — > LH) « 9 ~ 12.5 MeV depending on the choice of the LDMEs and 
T(hb — > LH) = 52.I3I728 keV. We conclude that for the bb system 0(a s v 2 ) corrections are 
not as important as in the cc system. We have also compared our theoretical results with 
experimental data || |j] and found that in general our calculations are consistent with data 
within theoretical and experimental uncertainties. In particular, we note that the operator 
evolution method used by us may provide a quite good estimate for the color-octet LDMEs. 
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Appendix A: EVOLUTION OF NRQCD MATRIX ELEMENTS C^sf 1 ) and ^(^J) 81 ) 
AT 0{a s v 2 ) 

In order to cancel the infrared divergence in short- distance coefficients corresponding to 
l P± Fock state, we should evaluate the NRQCD four-fermion operator Cf^S^ 8 ') and P^S^ 8 ') 
to sufficient orders. 

The 0(a s ) diagrams of LDMEs are sorted into three sets: self-energy diagrams which 
are related to self-energy corrections of external heavy (anti-) quarks, coulomb diagrams 
where the gluon is connected with both initial or final heavy quark and anti-quark and the 
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FIG. 9: The one-loop NRQCD diagrams which involve the Feynman rules up to 0(v 2 ). The 
Coulomb interactions and the cross diagrams have been suppressed. 

intersecting diagrams where the gluon is related to an initial heavy (anti-) quark and a final 
(anti-)quark. The results of the first two sets have been given in Refs. |I2] , pTjj and here we 
only calculate the intersecting diagrams which relate to the transition from 5* wave to P 
wave. 

Using the Lagrangian shown in Eq. @ and (f|) we can write the amplitudes of each 
diagram, the integrals of the diagrams shown in Fig. |9] are (other crossed diagrams are not 
shown) 



a f d D l q-q'-{q-l){q'-l)/l 2 1 - q 2 /2m 2 Q - q' 2 /2m* Q 

a+b+c Ws J (2ir)° mlill -l 2 + ie) k-I n -Ml + ie \ \a' — in — - 



2 

(Ala) 

/d = ^ 2 / ' J9^]d'i ; (q-i) 2 ; • ir / J ( 9 '-o 2 , • i ' ( Alb ) 

./ W [ go -z -i2^ + le ][^_j _i£_^ + ie ] 

where g = (go, Of) is the heavy quark external momentum and I = (Iq,1) loop integral 
momentum. Since there is no pole on the upper-half of Zn' s complex plane, the second 
integral yields zero. Contour- integrating the first integral over l around the Iq = \l\ — ie 
pole, we find 

d D -H q-q'-(q-l)(q>-l)/l 2 1 - q 2 /2m 2 Q - q' 2 /2m 2 Q 



1 a+b+c — 9 s 



(2*)"- 2m%\l\ Hi | _ ^ + £ + i£]Hi | __£_ + £!+ i£] 



(A2) 

Here we should expand the relative momentum in the denominator. Assuming q-l/rriQ, 
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q'-l/rriQ and I 2 /rriQ is far small than \l ][p8|, we get the required expansion 



9t f d D -H q-q'-(q-l)(q'-l)/l< 



2m 2 Q J (2n) D - 1 
x(l + (^) 2 + ^^ ^ 



(1 - q 2 /2m 2 Q - q' 2 /2m 2 Q ) 



\l\m Q \l\rriQ 
This integral can be reduced by taking the following substitution 

1 



+ (high order or irrelevant expansions). (A3) 



VI 3 -> 



D - 1 



-5 ij l 2 , 



VlH k V 



.(5 ij S kr + S ik 8 jr + S ir S kj )l 4 



(D-1)(D + 1) 

where S lJ is D — 1 dimensional Euclidean delta symbol. The integral yields 



2m^ 7T 2 D — 1 
Summing up all the diagrams we get 



D + 1 2mi 



1 



2a> } g - 2 
irm 2 ^ D — 1 



:i-§-|A(9 2 + g' 2 )) 

D + 1 2777,7} 



C(^ 



Recalling the definitions of O^pj 11 ) and P( 1 P 1 [1] ) we can write 

2(D-2)ai b) ( 1 1 



P|0( 1 4 8] )|i/) = (//|0( 1 ^)| J ff)Born + 



2A, 



(i/|P( 1 4 8l )|^) = (^|P(H 8| )|^)Born + 



2{D - 2)a 



(A4a) 
(A4b) 



(A5) 



(A6) 



(A7a) 



(D - l)7im 2 \e uv e IR 



(H\V{}PP))\H) 



C, 



2N r 



(A7b) 



Here we omit 0( x Pl B] ) and P^pf 1 ) coefficients since they are irrelevant in our case. The 
presence of UV divergence indicates that the LDMEs need renormalization. The relevant 
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counter-term in the MS scheme can be chosen 



H\OCS [ o ] )\H) = (H\OQsf)\H)^ + 



4a, 



3irmQ \euv 



+ ln47r - 7£ 



(H\Q(^ 1] )\H) 3 (H\VePn\H) 

^F 7 — — 



2A^ 



5m^ 



(H\veS®)\H) 



-2e 



{H\VCS [ o ] )\H) M + 



4a, 



37rmQ V e[/y 



2iV c 
1 



(A8a) 



+ ln47r - 'Je 



x C, 



(H\V(^l 1] )\H) 

2N r 



(A8b) 



where //a is the NRQCD renormalization scale. Combining Eq. ( |A7|) and ( |A8| ) we find 



(H\OCS [ * ] )\H) Bom = ^(H\OCS [ * ] )\H) 



(ma) 



+ 



4a, 



37rm^ V e IR 



h ln47r — 7^ 



x 



(^(^I^Born = V-A e (H\nX l )\Hr A> + 



-2e, 



2c 



(gjg(jf)jg) 3 (P^Pf 1 )^) 

Os ——z — k-L^f — 



2N r 



5m^ 



4a , 



3nm 2 \e IR 



2N r 



+ ln47r - 7 B 



(A9a) 



x [ ^ C F 



2c 



{H\VCP?)\H) 
2N r 



(A9b) 



Combining the self-energy part (see Eq. (B14) in Ref. [[[)) we get the total loop corrections 
of NRQCD LDMEs 



(H\OCS [ *)\H) Bom = ^(H\OCS [ * ] )\H)(^ + 



4a, 



1 



SimiQ \e IR 



2c 



(H\VCsf ] )\H) Born = ^(H\V(X l )\Hr A) + 7^ ( — + ^ ~ ^ 



A 6 2 -2 
4A C 

-2e, 



Of — — — K-Lyp — 



ln47T — 7s 



2N r . 



5m 2 Q 



2N r . 



(AlOa) 



(H\VCS l o ] ))\H) 



(ma) 



4a, 



3vrm^ \e IR 



2c 



2A, 



C, 



(AlOb) 
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These results are equivalent to Eq. (|3lD up to 0(a s v 2 ). 
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